output_version <- '2024-07-19'

## fishnet
fishnet_lake_chad_ntl_adm0_x_d01_shp <- paste0(wkdir,'/local/gis_data/003_boundaries/fishnet/fishnet_lake_chad_ntl_adm0_x_d01')
fishnet_adm0_only <- read_sf(dirname(fishnet_lake_chad_ntl_adm0_x_d01_shp),basename(fishnet_lake_chad_ntl_adm0_x_d01_shp))
st_crs(fishnet_adm0_only) <- "EPSG:4326"

##
## combined suitability FAO-FGGD 
##

# land suitability
landsuit_fn <- file.path(wkdir,'local/gis_data/007_environment/landsuitability/landsuitability.tif')
xfun::dir_create(dirname(landsuit_fn))

if(!file.exists(landsuit_fn)){
  url0 <- 'https://data.apps.fao.org/catalog/dataset/f8906e70-8922-11db-b9b2-000d939bc5d8/resource/2520786c-c57f-4d0c-8ccf-6bf1549d471f'
  url <- 'https://storage.googleapis.com/fao-maps-catalog-data/uuid/f8906e70-8922-11db-b9b2-000d939bc5d8/resources/Map6_62.zip'
  # Send the request
  response <- httr::GET(url, 
                        httr::write_disk(file.path(tmpdir, basename(url)),
                                         overwrite = TRUE))
  # Check the response status
  if (httr::status_code(response) == 200) {
    # Extract the file name from the URL
    downloaded_file <- file.path(tmpdir,basename(url))
    
    # Unzip the downloaded file to a temporary directory
    unzip(downloaded_file, exdir = tmpdir)
    
    # Example: List files in the temporary directory
    dir(tmpdir)
    
    cat("File downloaded successfully!")
  } else {
    cat("Error downloading file:", http_status(response)$message)
  }
  
  landsuit <- terra::rast(file.path(tmpdir,'comb_low'))
  writeRaster(landsuit,landsuit_fn)
} else {
  landsuit <- terra::rast(landsuit_fn)
}


is.factor(landsuit)
factor_raster <- landsuit

#Combined suitability of global land area for pasture and rainfed crops (low input level)
#Source: https://data.apps.fao.org/catalog/iso/f8906e70-8922-11db-b9b2-000d939bc5d8
#For each cell (same cell, 50km, 100km), we need the area shares of the various categories below:
  
# Specify your categories; you can customize this based on your needs
rat <- levels(factor_raster)[[1]]
dim(rat)
new_labels <- c('Inland Water',
                'unsuited: crops=no past=no',
                'marginal: crops=poor past=very poor',
                'marginal: crops=poor past=poor',
                'marginal: crops=poor past=suited',
                'good: crops=suited past=possible',
                'good: crops=wellsuited past=possible',
                'prime crops=prime past=possible'
)
#'Special Area'
rat$DESCRIPTION <- new_labels
levels(factor_raster) <- rat
levels(factor_raster)

require(exactextractr)

sum_cover <- function(x){
  list(x %>%
         group_by(value) %>%
         summarize(total_area = sum(coverage_area)) %>%
         mutate(sh_landsuit = total_area/sum(total_area)))
}


##
## RADIUS BUFFER 0 (grid cell only) 50KM (FROM POINT), 100KM (FROM POINT)
##

## memory issues
n_objectid_list=unique(fishnet_adm0_only$OBJECTID)
total_chunks <- split(n_objectid_list, ceiling(seq_along(n_objectid_list)/1000))
length(total_chunks)
library(doParallel)

for(chunk_index in 1:length(total_chunks)) {
  print(chunk_index)
  landsuit_chunk_out_fn <- file.path(wkdir,'proc_data','landsuit_chunks',sprintf('LANDSUIT_share_%s_of_%s_%s.csv',chunk_index,length(total_chunks),output_version))
  xfun::dir_create(dirname(landsuit_chunk_out_fn))
  
  if(!file.exists(landsuit_chunk_out_fn)){
    n_objectid_list_chunk <- as.numeric(unlist(total_chunks[chunk_index]))
    n_cores <- min(6,length(n_objectid_list_chunk))
    
    cl <- makeCluster(n_cores)
    registerDoParallel(cl)
    
    results <- foreach(objectid_j=n_objectid_list_chunk,
                   .export=c(ls()),
                   .packages = c("dplyr", "sf","exactextractr","spatialEco")) %dopar% {
                # Modulus operation
                if(objectid_j %% 100==0) {
                                     # Print on the screen message
                                     cat(paste0("iteration: ", objectid_j, "\n"))
                }
                
                fishnet_lake_chad_ntl_adm0_x_d01_shp_fn <- paste0(wkdir,'/local/gis_data/003_boundaries/fishnet/fishnet_lake_chad_ntl_adm0_x_d01.shp')
                fishnet_adm0_only_j = read_sf(fishnet_lake_chad_ntl_adm0_x_d01_shp_fn, query = sprintf('SELECT * from fishnet_lake_chad_ntl_adm0_x_d01 WHERE OBJECTID = %s',objectid_j))
                st_crs(fishnet_adm0_only_j) <- "EPSG:4326"
                
                if(dim(fishnet_adm0_only_j)[1]>0){
                  
                  fishnet_adm0_only_buf100km_bbox <- round(st_bbox(spatialEco::geo.buffer(x=st_centroid(fishnet_adm0_only_j), r=150000)),1)
                  
                  ## crop raster to 100km buffer
                  factor_raster <-  terra::rast(landsuit_fn)
                  factor_raster_crop <- terra::crop(factor_raster,fishnet_adm0_only_buf100km_bbox,snap="out")
                  
                  ## 0km
                  x <- exact_extract(factor_raster_crop, fishnet_adm0_only_j, coverage_area = TRUE, summarize_df = TRUE, fun = sum_cover)
                  
                  #merge the list elements into a df
                  landsuit_share0 <- bind_rows(x, .id = "objectid") %>%
                    mutate(buffer_km = 0)
                  
                  rm(x)
                  # projected 50km
                  fishnet_adm0_only_buf50km <- spatialEco::geo.buffer(x=st_centroid(fishnet_adm0_only_j), r=50000)
                  x <- exact_extract(factor_raster_crop, fishnet_adm0_only_buf50km, coverage_area = TRUE, summarize_df = TRUE, fun = sum_cover)
                  
                  #merge the list elements into a df
                  landsuit_share50 <- bind_rows(x, .id = "objectid") %>%
                    mutate(buffer_km = 50)
                  
                  # projected 50km
                  fishnet_adm0_only_buf100km <- spatialEco::geo.buffer(x=st_centroid(fishnet_adm0_only_j), r=100000)
                  x <- exact_extract(factor_raster_crop, fishnet_adm0_only_buf100km, coverage_area = TRUE, summarize_df = TRUE, fun = sum_cover)
                  
                  #merge the list elements into a df
                  landsuit_share100 <- bind_rows(x, .id = "objectid") %>%
                    mutate(buffer_km = 100)
                  
                  out <- rbind(landsuit_share0,landsuit_share50,landsuit_share100)
                  
                  # add objectid
                  out$objectid <- objectid_j
                  
                  # objectid id by iteration
                  out <- out %>%
                    mutate(objectid = objectid_j) %>%
                    dplyr::select(-total_area)
                  
                  # clean-up
                  rm(x,fishnet_adm0_only_buf100km,fishnet_adm0_only_buf50km)
                  
                  return(out)  
                }

                }                
    parallel::stopCluster(cl)
    
    system.time(landsuit_dt <- data.table::rbindlist(results))
    landsuit_df <- landsuit_dt %>% as.data.frame()
    data.table::fwrite(landsuit_df,landsuit_chunk_out_fn)
    
    # clean-up
    rm(landsuit_df,landsuit_dt)
    }
  
  
}

# merge all chuncks
landsuit_df_files <- sort(list.files(dirname(landsuit_chunk_out_fn),full.names=T))
tables <- lapply(landsuit_df_files, function(z) read.csv(z))
landsuit_df_all <- data.table::rbindlist(tables)

# change name for 
names(rat)[2] <- "landsuitcat"

landsuit_out <- merge(landsuit_df_all,rat,by.x='value',by.y='VALUE',all.x=T) %>%
  dplyr::select(objectid,sh_landsuit,landsuitcat,buffer_km) %>%
  arrange(objectid,landsuitcat,buffer_km)

landsuit_out_fn <- file.path(wkdir,'proc_data',sprintf('LANDSUIT_share%s.dta',output_version))

# Adding the label value
attr(landsuit_out$objectid, "label") <- "objectid of grid"
attr(landsuit_out$sh_landsuit, "label") <- 'Share of area covered by Land suitability category : FAO-FGGD'
attr(landsuit_out$landsuitcat, "label") <- 'Land suitability category : FAO-FGGD '
attr(landsuit_out$buffer_km, "label") <- 'Buffer'
attr(landsuit_out, "label") <- sprintf("sh_area_landsuit_out.R last run on %s",Sys.Date())

## Save the data
haven::write_dta(landsuit_out, landsuit_out_fn)
